DC-transport in superconducting point contacts: a full counting statistics view 
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We present a comprehensive theoretical analysis of the dc transport properties of superconduct- 
ing point contacts. We determine the full counting statistics for these junctions, which allows us 
to calculate not only the current or the noise, but all the cumulants of the current distribution. 
We show how the knowledge of the statistics of charge transfer provides an unprecedented level of 
understanding of the different transport properties for a great variety of situations. We illustrate 
our results with the analysis of junctions between BCS superconductors, contacts between supercon- 
ductors with pair-breaking mechanisms and short diffusive bridges. We also discuss the temperature 
dependence of the different cumulants and show the differences with normal contacts. 

PACS numbers: 74.50.+r, 72.70.+m, 73.23.-b 



I. INTRODUCTION 

The current-voltage (I-V) characteristics of supercon- 
ducting contacts have been the subject of investigation 
during the last four decades. The first experimental anal- 
yses were performed in tunnel junctions jlj. In this case 
the current inside the superconducting gap is suppressed, 
and the results can be accurately described with the BCS 
theory ,2]. However, very often a significant current is ob- 
served in the subgap region, which cannot be explained 
with the simple tunnel theory. The first anomalies were 
reported by Taylor and Burstein |3|] who noticed a small 
onset in the current when the applied voltage V was 
equal to the energy gap, A/e, in a tunneling experiment 
between two equal superconductors. Relatively soon af- 
terwards it was apparent 0, 01 that not only is there 
an anomaly in the current at eV = A, but in fact at 
all submultiples 2A/n, where n is an integer. This set 
of anomalies is referred to as subharmonic gap structure 
(SGS), and its temperature and magnetic field depen- 
dence were thoroughly characterized [a, 13 • 

The first theoretical attempt to explain the SGS was 
done by Schrieffer and Wilkins Q, who noticed that if 
two electrons could tunnel simultaneously, this process 
would become energetically possible at eV = A, and 
cause the structure in the I-V observed by Taylor and 
Burstein [2|. Within this multiparticle tunneling theory 
the origin of the SGS would be the occurrence of multiple 
processes in which n quasiparticles cross simultaneously 
the contact barrier. The original perturbative analysis of 
this theory has serious problems. In particular, the cur- 
rent was found to diverge at certain voltage, which avoids 
to calculate meaningful I-Vs within this approach. A 
second explanation was put forward by Werthamer [lfj, 
who suggested that the SGS could be caused by a self- 
detection of the ac Josephson effect. The main problem 
of this explanation is that it invokes two different mecha- 
nisms for the odd and even terms, while the experimental 
current jumps are identical for both series. In 1982 Klap- 



wijk, Blonder and Tinkham [ll| introduced the concept 
of multiple Andreev reflection (MAR) . In this process a 
quasiparticle undergoes a cascade of Andreev reflections 
in the contact interface (see Fig. QJ. They showed that 
a MAR in which a quasiparticle crosses the interface n 
times becomes possible at a voltage eV — 2A/n, which 
explains naturally the SGS. The quantitative analysis of 
the I-Vs was based on a semiclassical app roach which 
fails away from perfect transparency |12lll3l |. A few years 
later, Arnold reported the first fully microscopic calcula- 
tion of I-Vs based on a Green's function approach |l4| . 

The theoretical discussion was finally clarified with the 
advent of modern mesoscopic theories. Using the scatter- 
ing formalism [JJJ [l|| [jj| and the so-called Hamiltonian 
approach [||| , different authors reported a complete anal- 
ysis of the dc and ac Josepshon effect in point contacts. 
These theories clearly showed that the MARs are respon- 
sible of the subgap transport in these systems. They also 
showed that the multiparticle tunneling of Schrieffer and 
Wilkins and the MARs are indeed the same mechanism. 
The new microscopic theories have also allowed the cal- 
culation of a series of properties such as resonant tunnel- 
ing [jjj|2(|, shot noise |2lll22| and the Shapiro steps |23| . 

From the experimental point of view, the main prob- 
lem has always been the proper characterization of the 
interface of the superconducting contact. Uncertainties 
in the interfaces properties often avoid a proper com- 
parison between theory and experiment. The situation 
has considerably improved with the appearance of the 
metallic atomic-sized contacts, which can be produced 
by means of scanning tunne ling microscope and break- 
junction techniques [H |H M, IH |H El H3, IH H2 . 
These nanowires have turned out to be ideal systems to 
test the modern transport theories in mesoscopic super- 
conductors. Thus, for instance Scheer and coworkers J23 
found a quantitative agreement between the measure- 
ments of the current-voltage characteristics of different 
atomic contacts and the predictions of the theory for a 
single-channel superconducting contact [la. [l8j . These 



experiments not only helped to clarify the origin of the 
SGS, but also showed that the set of the transmission 
coefficients in an atomic-size contact is amenable to mea- 
surement. This possibility has recently allowed a set of 
experiments that confirm the theoretical pred ictions for 
transport properties such as supercurrent j^Jl, noise [32 ] 
and even resonant tunneling in the context of carbon 
nanotubes 33]. From these combined theoretical and 
experimental efforts a coherent picture of transport in 
superconducting point contacts has emerged with multi- 
ple Andreev reflections as a central concept. 

The most recent development in the understanding of 
the dc transport in superconducting contacts is the anal- 
ysis of the full counting statistics |3J, [3!| . Full counting 
statistics (FCS) is a familiar concept in quantum optics 
(see for instance |36J]). which has been recently adapted 
to electron transport in mesoscopic conductors by Levi- 
tov and coworkers |37| . FCS gives the probability P(N) 
that iV charge carriers pass through a conductor in the 
measuring time. Once these probabilities are known one 
can easily compute not only the mean current and noise, 
but all the cumulants of the current distribution. Since 
the introduction of FCS for electronic systems, the the- 
ory has been sophisticated and applied to many different 
contexts (see Ref. 38] for a recent review). 

The counting statistics of a one-channel quantum con- 
tact has the surprisingly simple form of a binomial dis- 
tribution P(N) = (%)T N (1 - T) M - N , where T is the 
transmission pro bability and M ~ V is the number of 
attempts [33, [33 • The generalization to many contacts 
and/or finite temperatures is straightforward, by noting 
that different energies and channels have to be added 
independently. In this way, the counting statistics of 
diffusive contacts at zero temperature [40] and at finite 
temperatures |4l| could be obtained usin g th e universal 
distribution of transmission eigenvalues |42l |43| . It is 
worth to note, that the FCS in the limit of small trans- 
parency reduces to a Poisson distribution, which can also 
be obtained using classical arguments and neglecting cor- 
relations between the different transfer events. Interest- 
ingly, the Poissonian character allows to directly extract 
the charge of the elementary event, which can be used to 
determine e.g. fractional charges J44|,[45|, 46j . A general 
approach to obtain the counting statistics of mesoscopic 
condutors was formulated by Nazarov J4l| using an ex- 
tension of the Keldysh-Green's function method, which 
allowed to present the counting statistics of a large class 
of quantum contacts in a unified manner [47| • In Ref. J34J 
we have shown, how this method can be used for a time- 
dependent transport problem like a superconducting con- 
tact out of equilibrium. 

The counting statistics of a contact between a nor- 
mal metal and a superconductor at zero temperature 
and eV <C A was shown to be again binomial with the 
important difference that only even numbers of charges 
are transferred [48|. The probability of an elementary 
event is then given by the Andreev reflection coefficient 
R A = T 2 /(2-T) 2 H<|. Again, the generalization of this 



result to many channel conductors is obtained by sum- 
ming over independent channels. For a diffusive metal 
the resulting statistics was shown to be an exact replica 
of the corresponding statistics for normal diffusive trans- 
port, provided the double charge transfer is taken into ac- 
count |50|. This holds for coherent transport eV <C £r/i, 
where Eth is the inverse diffusion time, as well as in the 
fully incoherent regime eV 3> Etk |51| . For intermedi- 
ate voltages, correlations of transmission eigenvalues at 
different energies modify the distribution of transmission 
eigenvalues 521, which lead to a nonuniversal behavior of 
the transport statistics, predicted theoretically [53] and 
confirmed experimentally [54j . Here, we note that a dou- 
bling of the noise was experimentally observed in diffusive 
wires [55] , confirming earlier theoretical predictions [5 6) . 
However, to trace this back to a doubling of the elemen- 
tary charge transfer follows only from an analysis of the 
counting statistics. A direct experimental determination 
of the doubled charge transfer was recentl y pe rformed in 
a conductor containing a tunnel junction [57[. Here, the 
underlying statistics is Poissonian and the noise dir ectly 
gives access to the charge of the elementary event |5a.l59| . 

An interesting problem occurs, when one applies the 
concept of counting statistics to a supercurrent through a 
quantum contact [47j • The resulting statistics can not be 
directly related to a probability distribution, since some 
of the 'probabilities' would be negative. A closer inspec- 
tion of the formalism showed, that the interpretation of 
probabilities relies on the prop er definition of a quantum 
measuring device |60t l(3ll l62j| . As we will see below, in 
superconducting contacts out of equilibrium these prob- 
lems do not occur and all probabilities are positive. 

In Ref. J34| we have demonstrated that the charge 
transport in superconducting point contacts out of equi- 
librium can be described by a multinomial distribution of 
processes in which a multiple charge is transferred. More 
importantly, we have shown that the calculation of the 
FCS allows us to identify the probability of the individ- 
ual MARs and the charge transferred in these processes. 
This information probably provides the deepest insight 
into the transport properties of these systems. In this 
sense, in this work we present a comprehensive analysis 
of the dc transport properties of superconducting point 
contacts from the point of view of the FCS. We shall show 
that even in the most well-studied situations the FCS 
provides a fresh view. Moreover, we show that the FCS 
allows a unified description of many different type con- 
tacts. We also extend the analysis presented in Ref. J34J 
to finite temperature. 

The paper is organized as follows. In section II, af- 
ter introducing some basic concepts of charge statis- 
tics, we discuss the calculation of the cumulant gener- 
ating functional within the Keldysh-Green's function ap- 
proach. Section III is devoted to the calculation of the 
MAR probabilities at zero temperature. We present both 
the results of a toy model and the full expressions. In sec- 
tion IV we apply the results of the previous section to de- 
scribe the different transport properties of three different 




logarithm of the characteristic function and is defined by 



eV> 2A/4 



FIG. 1: Schematic representation of the MARs for BCS su- 
perconductors with gap A. We have sketched the density of 
states of both electrodes. In the upper left panel we describe 
the process in which a single electron tunnels through the sys- 
tem overcoming the gap due to a voltage eV > 2 A. The other 
panels show MARs of order n — 2,3, 4. In these processes an 
incoming electron at energy E undergoes at least n — 1 An- 
dreev reflections to finally reach an empty state at energy 
E + neV. In these MARs a charge ne is transferred with a 
probability, which for low transparencies goes as T n . At zero 
temperature they have a threshold voltage eV — 2A/n. The 
arrows pointing to the left in the energy trajectories indicate 
that a quasiparticle can be normal reflected. The lines at en- 
ergies below E and above E + neV indicate that after a detour 
a quasiparticle can be backscattered to finally contribute to 
the MAR of order n. 



situations: (i) a contact between BCS superconductors, 
(ii) a contact between superconductor with a modified 
density of states due to a pair-breaking mechanisms, and 
(iii) a short diffusive SNS contact. In section V we an- 
alyze the transport at finite temperature paying special 
attention to the third cumulant. Finally, we present our 
conclusions in section VI. 



II. DESCRIPTION OF THE FORMALISM 
A. Some basic concepts 

Our goal is to calculate the full counting statistics of 
a superconducting contact. This means that the quan- 
tity that we are interesting in is the probability P to (N), 
that N charges are transferred through the contact in the 
time interval to- Equivalently, we can find the cumulant 
generating function (CGF) St (x) , which is simply the 



exp(5 t0 (x))=J2 P *o (N) exp( ? iV x ) 



(1) 



A' 



Here, x is t ne so-called counting field. From the knowl- 
edge of the CGF one easily obtains the different cumu- 
lants that characterize the probability distribution 



Qn 



(2) 



x=o 



Notice that the first cumulants are related to the mo- 
ments of the distribution as follows 



Ci = N = J2NP t0 (N) , C 2 = (N-Nf, 



N 



C 3 = (N- N) 3 , d = (N - N) 4 - 3C£ , (3) 

and so on. It is also important to remark that these cu- 
mulants have a simple relation with the relevant trans- 
port properties that are actually measured. Thus, for 
instance, the mean current is given by / = (e/to)Ci 
and the symmetrized zero frequency noise is given by 
Si = (2e 2 /to)C2UM- For higher cumulants such rela- 
tions are not straightforwardly obtained, but it can be 
shown that the cumulants defined above correspond to 
the observable quantities in an electron counting exper- 
iment J43, Lm l&M- Thus, the cumulants represent all 
information, which is available in a measurement of the 
charge accumulated during the observation period to- 



B. Keldysh-Green's function approach to FCS 

As mentioned above, our system of interest is a voltage- 
biased superconducting point contact, i.e. two supercon- 
ducting electrodes linked by a constriction, which is much 
shorter than the superconducting coherence length. We 
concentrate ourselves on the case of a single channel con- 
tact described by a transmission probability T . The main 
difficulty in the determination of the FCS arises from the 
ac-Josephson effect. Here, a constant applied bias volt- 
age eV gives rise to time-dependent currents as a conse- 
quence of the Josephson relation {d/dt)<p(t) — 2eV/h. In 
the long-time limit to ^> h/eV these oscillating currents 
do not contribute to the net charge transfer, in which we 
are interested. However, this intrinsic time-dependence 
is reflected in the CGF and a little care has to be taken, 
when the FCS is defined. 

To obtain the FCS in a superconducting point con- 
tact we make use of the Keldysh-Green's function ap- 
proach to FCS introduced by Nazarov and one of the 
authors |4l|, |42|, and we refer to reader to these pa- 
pers for further details on the basis of this theoretical 
approach. In what follows, we concentrate ourselves on 
the specific difficulties introduced in the case of a contact 
between two superconductors. Our starting point for the 



determination of the CGF is to define the relation be- 
tween the CGF and the counting current in analogy to 
Refs. EH!!: 



JU (x) = -/ <im\.n. 



(4) 



This scalar current can be calculated in terms of the ma- 
trix current which describes the transport properties of 
the contacts. Nazarov has shown that, in the case of 
short junctions the matrix current (in Keldysh-Nambu 
space) adopts the following form [22J 



2T[G 1 ( X )?G 2 ] 



(t,0 



i(x,t,t') = — 

7T \ v 4 + T({G 1 ( X )?G 2 }-2) / 

(5) 
Here Gi^ 2 )(t,t') denote the matrix Green's functions on 

the left and the right of the contact. In our problem 
these functions depend on two time arguments and the 
products ® appearing in Eq. (JSJ) should be understood as 
convolutions over the intermediate time arguments, i.e. 
(A®B)(t,1?) = Jdt"A(t,t")B(t",t'). It is worthwhile 
to n ote, that the derivation for the matrix current in 
Ref. [&| was done for Green's functions in the static situ- 
ation, in which case all Green's functions depend only on 
t — t' . However, the derivation can be directly taken over 
to time-dependent problems, because the time-dependent 
Green's functions satisfy the normalization condition 



{G®G)(t,t') = S{t-t') 



(6) 



Finally, the time-dependent scalar current is obtained 
from the matrix current by 



I( X ,t) = ^Tr[f K I(x,t,t)] , 



(7) 



where tk = &3T3 is a matrix in Keldysh(")-Nambu(~) 
space. <7i(Tj) are the standard Pauli matrices in 
Keldysh(Nambu)-space. 

Let us now describe Green's functions entering Eq. J5J . 
The counting field x 1S incorporated into the matrix 
Green's function of the left electrode as follows 

Gi (x, t, t') = e- lxfK / 2 G x (t, t')e lxfK/2 . (8) 

Here Gi (£,£') is the reservoir Green's function in the ab- 
sence of the counting field. We set the chemical potential 
of the right electrode to zero and represent the Green's 
functions by 

Gx(M') = e^ (t)?3/2 Gs(t-t')e" #( *'^ 3/2 (9) 

and G 2 (M') = Gs(£ - £')■ Here : 0(*) = 4>o + (2eV/h)t 
is the time-dependent superconducting phase difference, 
and <po is its dc part. Gs is the Green's function of a 



superconducting reservoir (we consider the case of a sym- 
metric junction), which reads 



G s (t - 

G S (E) 



iE(t-t') 



dE G s {E)e 



(A -R)f + R (A- R)f 
(A-R)(l-f) (R-A)f + A 



(10) 



Here, R(A)(E) are retarded and advanced Green's func- 
tions of the leads and f(E) is the Fermi function. Ad- 
vanced and retarded functions in 11U|) have the Nambu- 
structure R(A) — g R ' A T3 + / R,A fi fulfilling the normal- 
ization condition f 2 + g 2 = 1 . They depend on energy 
and the superconducting order parameter A. 

Using the time dependence of the leads Green's func- 
tions it is easy to show from Eq. J5J) that the scalar cur- 
rent admits the following Fourier series 



I( X ,t)=^I n (x)e 



iri(p(t) 



(ii) 



which means that the current oscillates with all the har- 
monics of the Josephson frequency. It is important to 
stress that the components I n (x) are independent of dc 
part of the superconducting phase. In this work we only 
want to consider the dc part of the CGF. For this pur- 
pose, we take the limit of a long measuring time to, much 
larger than the inverse of the Josepshon frequency, and 
hereafter we drop the subindex to in the expression of 
the CGF. From Eq. @ and Eq. JT1J it is obvious that 
by selecting the dc component, the dc part of the phase 
drops the calculation and the CGF is free of the problems 
related to g au g e invariance found for the dc Josephson ef- 
fect miail. 

Keeping in mind the presence of the time integration 
described aboved, and with the help of Eqs. I|5I7|) . one 
can integrate Eq. J3J to obtain the following expression 
for the CGF of superconducting constrictions [43 



S(x) = -j-Tr In 



l + j({Gi(x),G 2 } (3 -2) 



(12) 



The symbol ® implies that the products of the Green's 
functions are convolutions over the internal energy argu- 
ments, i. e. 

(G 1 ®G 2 )(E,E')= I dE 1 G 1 (E,E 1 )G 2 (E 1 ,E r ). (13) 

The trace runs not only over the Keldysh-Nambu space, 
but also includes a trace in the energy arguments, i.e. 
JdE g(E,E). 

The time-dependent Green's functions of Eq. JSJ fulfill 
the normalization condition of Eq. ©. This enables us 
to use the relation 



2 - {Gi, G 2 } 8 — (Gi - G 2 ) 
to write the CGF as 

S( X ) = ^Tr{lnQ + +lnQ_} 



(14) 



(15) 



where Q± = 1 ± (VT/2) (Gi(x) - G 2 ). One can show 
that both logarithms give the same contribution, and 
therefore we concentrate in the analysis of the first one, 
and we drop the subindex +. Additionally, we use the 
relation Tr In Q = In dct Q to write the CGF as 



S(X) 



lndetQ(x). 



(16) 



Thus, at this stage the calculation reduces to the cal- 
culation of the determinant of a infinite matrix. Due to 
the time dependence of the lead Green's functions their 
form in energy space is G(E,E') = J2 n Go,n(E)6(E — 
E' + neV), where n = 0,±2. This implies that the ma- 
trix Q also admits the same type of representation, which 
in practice means that Q is a block-tridiagonal matrix of 



the form 



V 



Q 








-2,- 


-4 Q-2,-2 Q-2,0 





Q-2,0 Qo,a Qo,2 




Q2,o Q2.2 



J 



Q(E + neV,E- 



where we have used the notation Q rh 
meV). The different (4 x 4) matrices Q n ,m have the fol- 
lowing explicit form in terms of the advanced and re- 
tarded Green's functions g R ' A and f R ' A (remember that 
we consider a symmetric junction) 



J 



Q, 



( Pn+l - Pn + 9n±l ~ 9n 
' -Pn-fi 

e tx S n+ i - 5 n 

-S n 



-Pn ~ /* 



Pn ~ Pn-ig* ~ gn-1 



-e lx 5 n -i + S n 



e lX p n +l - Pr, 



Pn - Pn+l + 9n+l - 9n 



in,n+2 



-e lx p n -i + p r , 

-In + Pn 



/0 e-HPn+l + tf+l) pn+l \ 





e*(/£.i - Pn+l) 





Pn-l 







la+l 

\0 

/ 0\ 

e lX {~Pn~l +f,tl) Pn-l 



V S n -i e-^ifti - Pn-i) ) 



~fn + Pn Pn-l ~ Pn + 5„ ~ 9n-l I 



(17) 



r 



where we have used the shorthand notation g£ 



r.a _ 



,R,A 



(E + neV), p — (g A — g R )f, f being the Fermi 



function, p = (f A - /«)/, 6 

s = (f A -f R )(i-f). 



(g A - g R )(l - f), and 



One can restrict the fundamental energy interval to 
E — E' e [0, eV] , and therefore the CGF adopts the form 

S(x) = (W l )J eV dE lndet( 3- From E q- (113, it is ob- 
vious that det Q can be written as the following Fourier 
series in \ 



form 



Six) 



where 



eV 



dE \n 



1+ J2 PniE,V) ie in * - 1) 



n— — 00 



P n iE,V) = P n iE,V)/ Y, KiE,V). 



(19) 



(20) 



detQ(x) = 



n— 00 



P' n iE,V)e mx 



(18) 



where the coefficients P^iE^V) have still to be deter- 
mined. Keeping in mind the normalization S'(O) = 0, it 
is clear that one can rewrite the CGF in the following 



Eq. IJ19JI has the form of the CGF of a multinomial dis- 
tribution in energy space (provided more than one P n 
is different from zero). The different terms in the sum 
in Eq. I|19|) correspond to transfers of multiple charge 
quanta ne at energy E with the probability P n {E,V), 
which can be seen by the (27r/n)-periodicity of the ac- 
companying x-dependent counting factor. This is the 
main result of our work and it proves, that the charges 
are indeed transferred in large quanta. Of course, we still 
have to determine the probabilities P n iE, V), which is a 
non-trivial task and it will the goal of the next section. 



C. Cumulants 

As explained before, from the CGF one can easily cal- 
culate the cumulants of the distribution and in turn many 
transport properties. Of special interest are the first 



J 



three cumulants C\ , C2 and C3 , which correspond to the 
average, width and skewness of the distribution of trans- 
mitted charge, respectively. From Eq. J2J an d Eq. (TSJ, 
it follows that these cumulants can be expressed in terms 
of the probabilities P n (E, V) as follows 



Ci = 



C 2 = 



C* = 



f [ V dE I E n3p - + 2 (e nP >) - 3 (e nP «) 




(21) 
(22) 

(23) 



These expressions are a simple consequence of the fact 
that the charge transfer distribution is multinomial in 
energy space. At zero temperature the sums over n are 
restricted to positive values (n > 1). We remind the 
reader that the first two cumulants are simply related to 
the dc current, / = (e/tQ)Ci, and to the zero- frequency 
noise S/ = (2e 2 /t )C 2 - 

It is instructive to discuss some consequences of these 
expressions. Let us first recall, what happens when only 
one process contributes, which has, e.g., the order n. The 
first three cumulants are 



Cl-r 



cw 



eV 



t dE 



n h 

V 



p„ 



C 3; „ = n 3 f 
Jo 



t dE 



t dE 



Pn (1 - Pn) , 



(24) 
(25) 



P n (1 - P n ) (1 - 2P n ) . (26) 



We see, that the i th cumulant is proportional n\ i. e. 
the i th power of the charge of the respective elementary 
event. The expressions under the integral in Eqs. (|24I26|) 
have the same form as for binomial statistics, however in 



r 



general the P n {E, V) depend on energy in a nontrivial 
way and the energy-integrated expressions for the cumu- 
lants do not correspond to binomial statistics. A simple 
interpretation in terms of an effective charge transferred 
is only possible if P n (E, V)cl for all energies, in which 
case one recovers the standard result for Poisson statis- 
tics, Ci- n — 7i l_1 Ci ;n . According to Eq. (|2lfl) the sign of 
the spectral third cumulant can be positive or negative, 
depending on the size of P n (positive for P n < 1/2 and 
negative for P n > 1/2). The overall sign depends on the 
energy average and is not simple to predict. Note, how- 
ever, that the probabilities of MAR-processes of higher 
orders decrease approximately as T n . We may there- 
fore speculate that to obtain a negative third cumulant 
for higher order processes we will need more open con- 
tacts (a rough estimate is thus that T > 1/ \/2 to have 
P n > 1/2 and, therefore, C 3 < 0). 

The general statistics (|19H is a multinomial distribu- 
tion and it is therefore interesting to compare with inde- 
pendent binomial distributions. This is most easily done 
by assuming, that only two processes compete. Taking 
these processes to be of order n and m the first three 
cumulants read 



Cl\nm — C^n + C2:m ~ 2nm 



Cz\nm = C$- n + C^ m — 3nm 



eV 



t dE 



t dE 



P P 



P n P m [n(l - P n ) + m(l - P„ 



(27) 
(28) 

(29) 



r 



We see that the first cumulant is just the sum of the therefore have to look at higher cumulants to gain infor- 
contributions of the different processes n and m and we 



mation on correlations between the processes of different 
order. In both, the second and the third cumulant, such 
correlations appear and it is evident from Eqs. (|28|1 and 
(|29[l that both are reduced below the value obtained for 
independent binomials. The correlation terms appear in- 
side the energy integration and therefore both processes 
must be possible at the same energy. 

Finally, we note that in order to study correlation be- 
tween N different processes one would have to look at the 
Nth order cumulant. This becomes clear if one notices 
that only the TVth cumulant contains a term with prod- 
ucts of N probabilities and therefore the possibility to 
have a product of probabilities of N different processes. 



super: 


E> A 


|A| >E 


-A > E 


normal: 


E>eV 




E <eV 


5n (x) 


K-ix) 





K+{X) 


322 (x) 


-K-(-x) 





-K+(-x) 


912(21) (X) 





e ±lXT3 






TABLE I: Green's functions in the toy model. The indices 
g a p denote the respective element in Nambu space. K± = 
=Fr3 — 2f±e x denotes a matrix in Keldysh space. The table 
holds for left and right terminal, provided the energies and 
the counting fields are chosen properly. 
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III. MAR PROBABILITIES: 
TEMPERATURE 



ZERO 



This section is devoted to the calculation of the proba- 
bilities P n (E, V) at zero temperature. First, we discuss a 
simple model which nicely illustrates the transmission de- 
pendence of these probabilities, and secondly we present 
the general expressions. 



A. Toy model 

To obtain a feeling for the forthcoming calculations we 
will now study a strongly simplified model of a supercon- 
ducting contact. For that purpose, let us assume that 
we can neglect the Andreev reflections for energies out- 
side the gap region and replace the quasiparticle density 
of states by a constant for \E\ > A. Furthermore, we 
neglect that energy-dependent phase shift ~ acos(_E/A), 
usually associated with the finite penetration of excita- 
tions close to the gap edge. Mathematically, this means 
that we set f R > A {\E\ < A) = 1, g R{A \\E\ > A) = ±1, 
and both are equal to zero otherwise. This simplifies the 
calculation a lot, since the matrix Q in Eq. Hlbfl now 
becomes finite. In particular, for subharmonic voltage 
eV = 2A/n the matrix is also energy-independent. It 
is interesting to note that the toy model is also able to 
describe the counting statistics of normal contacts and 
Andreev contacts. 

To facilitate the discussion of the matrix structure it is 
useful to introduce the 2cg)2 matrix in the Keldysh sub- 
space 



K±{x) = TT 3 -2f ± e 



±iX 



(30) 



where f^ are Pauli matrices and f± = (fi ± if 2) /2. In 
fact, K± correspond to occupied (empty) quasiparticle 
states (for E > |A|). The matrix structure for supercon- 
ducting or normal terminals is summarized in Table [I] 
The counting statistics is obtained from the general re- 



S(x,V) 



2io 

h 



dE Tr In 



l + : hr{G l { x )-G 2 ) 



(31) 

The incorporation of the energy discretization is obtained 
by a redefinition of the trace in the above formulas, and 
a limitation of the energy integration to an interval of 
width eV. Note, that we have to evaluate only one of the 
two terms Q±, since the FCS can only depend on T and 
not on VT. 



To calculate the determinant we note that Q is a band 
matrix of width 3 in the energy index. Then the following 
reduction formula for the determinant is useful (assuming 
a block starts at some n, which we arbitrarily set to zero): 



Qo,Q Qo,2 

Q2 : Ql,2 Q2,4 



Q.4,2 Qaa 

'•■ '■ 

?2,2 — Q2,oQo,oQo,2 Q2,4 



(32) 



Qo,o 



t/4,2 




14,4 



Another useful property (which holds in the toy model) 
is the Nambu structure of the Q's, see Eq. (|17|) and Ta- 
ble HJ diagonal components in energy space, i. e. Q n ,n, 
are always block-diagonal in Nambu space and the off- 
diagonal components Q n .n±2 are purely off-diagonal in 
Nambu space and diagonal in Keldysh-space. Conse- 
quently, Qn-2,n-2 - Qn~2,nQn)nQn,n-2 appearing in tllC 

expansion of the determinant is block-diagonal again and 
the whole calculation of the determinant l|16l) boils down 
to a recursive calculation of determinants and inversions 
of 2x2-matrices. This will become more clear, when we 
will treat the explicit examples below. 
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1. Normal Contact 

It is instructive to demonstrate the procedure first for 
a normal contact. The Green's functions are (we restrict 
the calculation here to electron block, the hole block gives 
actually the same contribution) 



Ch 



K_( x ) , n>0 
K+( X ) , n<0 



Ch 



K-(0) , n>0 
K+(0) , n < 



,( 33 ) 

Note that we have chosen the fundamental energy inter- 
val [— eV/2, eV/2], since then the Green's functions are 
constant inside each interval. Then we find 



(Q-l)n, 

y/T/2 



= 0~n,r. 



f + (/*-l) 


n > 




T3 + t+e tx — f- 


n = Q . 


(34) 


f_(e- l x-l) 


n < 





The matrix Q has thus block diagonal form. The blocks 
n > and n < are tridiagonal and the determinants are 
all equal to 1. The remaining determinant of the n = 
block is 



det 



i + Vt y/re 



T 1 



= 1 - T + Te l 



(35) 



The CGF is finally S{\) = (2eVt a /h)\Ti(l + T(e^-l)) in 
agreement with Levitov and Lesovik 37]. Notice that a 
factor of 2 enters the CGF, because we get an additional 
contribution from the hole block (thus it is due to spin). 



2. Andreev contact 

We now consider a contact in which one of the sides is 
superconducting and the other is a normal metal. Again, 
the calculation can be done in a similar way. Here we 
apply a voltage \eV\ < A to the normal contact. The 
Green's functions are again diagonal in the energy space, 
since we assume that the superconductor is at zero poten- 
tial. For the normal metal we find (taking as fundamental 
energy interval [— eV, eV]) 



Ch 



Ch 



K-(x) , n>{) 

k+( x ) , n<0 ' 

-k_( X ) , n>0 

-K+(x) , n<0 



(36) 



and for the superconductor I G2 ) = I G2 ) = f 1 and 

otherwise. The only non-zero block is the n — energy 
block 



Gi(x)-G 2 



K-{X) -1 
-1 -K + (- X )J ' 



(37) 



which yields for the CGF in the form (|12|) the determi- 
nant of 



Q — \ 1 



1-2: 

± 2 



i(K+-Kh 



(K, 



K_) 



(38) 



To calculate the determinant we subtract from rows 3 and 
4 the rows I and 2 multiplied with ?(1 - ?)(■&+ - KJ) 



The matrix is then tridiagonal and its determinant is 



and make use of the fact that (K_ — K^ 



4(1- 



T 

1 

2 



1 



rpj, 



(2-T) 2 



,*2 X 



1) 



(39) 



The prefactor is canceled because we are operating un- 
der the In and have to normalize. Notice that the evalua- 
tion of the determinant outside the transport window can 
be done in a similar way. One obtains for the determi- 
nant of one block (1 - T/2) 2 -T 2 (K_(x) - k-(- X )) 2 = 
(1 — T/2) 2 , which is independent of the counting field x 
and is therefore canceled after normalization of the CGF. 
Finally we obtain for the FCS (collecting all prefactors) 

M 



s(x) 



2eVtn 



In 



f 



Tf 



(< 



i2 X 



1 



(40) 



The statistics corresponds to a binomial distribution of 
charge transfers. The Andreev reflection leads to a ir- 
periodicity in x which shows that only couples of charges 
can be transferred and the charge transfer probability for 
odd charge numbers vanishes. The number of attempts, 
determined by the prefactor of the In in (|40(l . remains 
unchanged in comparison to the normal case. 



3. Superconducting point contact 

We now come to the main subject of the article, a 
point contact between two superconducting banks held 
at different chemical potentials. To write down the gen- 
eral matrix structure of the FCS in the toy model, let us 
first obtain the condition for energies to be subgap. Here, 
we restrict ourselves to subharmonic voltages, which we 
write in general as eV = 2A/(N — 1), where N denotes 
the order. The dominating charge transport mechanism 
we expect is that N charges are transferred. In the toy- 
model, it is the only transport mechanism (since Andreev 
reflections above the gap are neglected). To obtain a 
single- valued matrix entries, it is favourable to choose as 
fundamental energy interval [0, eV] for even N = 2M 
and [-eV/2, eV/2] for odd N = 2M-1. For the Nambu 
row indices of the Green's function of the left terminal 
we find 



Nambu 


Order 


\E\ < A 


upper 


odd 


-M < n < M - 1 


lower 


odd 


-M + 1 < n < M 


upper 


even 


-M -Kn<M-l 


lower 


even 


-M <n< M 



(41) 



The row indices in Nambu space of the right Green's 
functions have the energy arguments of upper and lower 
row interchanged. 

To clarify the matrix structure we have prepared a 
small table. Each entry denotes the energy for the 



structure 



9u 






, where the second (Nambu-) in- 



dex i = l,2 plays no role. The entries are denoted by 
for E > A, for \E\ < A, and - for E < -A. 



n 


N = 2 


iV = 3 


A = 4 


A = 5 


N = 6 


2 


+ + 
+ + 


+ + 
+ + 


+ + 

+ + 


+ 
+ 


+ 
+ 


1 


+ + 
+ + 


+ 
+ 


+ 
+ 















+ 
+ 






















-1 


- 
- 


- 
- 

















-2 


- - 


- - 


- 
- 


- 
- 







-3 










- 




(42) 



We observe that the matrix structure in all cases is sim- 
ilar. A block with and + elements, i.e. connecting the 
quasiparticle states above the gap to the subgap region is 
followed a number of blocks inside the gap (depending on 
the applied voltage and, finally, is connected by a block 
with and — elements to quasiparticle states below the 
gap. 

Let us now discuss the case N — 2 (eV = 2A). Here 
the relevant 8 x 8-matrix is 



Q -l 

f/2 



/KM) 



\ 






-1 





K-(0) 

e Jxf3 




-1 

e~ lxf3 

-K+(0) 

-K+(-x) ) 



(43) 

We observe, that the matrix decouples into two blocks of 
4x4 matrices 



!2A 



= 1 



K-(X) 
-1 



-1 
-K + {-X) 



and 



Q2B-1 + —^ ei%f3 _ k+{Q) 



(44) 



(45) 



By comparison with Eq. I|37|l we see that lndetQ2A 
yields the counting statistics of usual Andreev reflec- 
tion. Qib gives actually the same result. This is most 
easily seen, if the unitary transformation IJQibU* with 



U = diag(e 



■T3X/2 



,e 



1T3 */ 2 ) is applied, which transforms 



Qib into Qia- Note, that the signs of the off-diagonal 
matrices do not matter, since they can be eliminated by 
similar unitary transformations. The counting statistics 
is therefore given by Eq. (|40|) . the same as for the An- 
dreev contact. 

Now we come to the slightly more complicated case 



N = 3 (eV = 2 A/2). Here we encounter the matrix 

10 \ 



( K-bt) o 

#_(0) e" lxf3 
e l * fa 

-10 



0-1 







-1 

e ixn -K+(0) 

-K+(-x) J 

Once again, the matrix decouples into two blocks (rows 
1,4,5 and rows 2,3,6). The first block is 




1 





e -iXT3 

K + (0) 



(47) 



It is already evident, that we will encounter a 
three particle process, if we apply the transformation 
U =diag(exp(ixf 3 ),exp(ixr 3 ),l). This yields 



UQ 3 aU j = 1 




-(3x) 

-1 


-1 





1 





1 


K + (0) 



(48) 



Evaluating the determinant we obtain the counting 
statistics (including the other block, see below) 



s(x) 



2eVt n 



In 1 



T^3 



(4 - 3T) 2 



(e-* 3 * - 1) 



(49) 



Evidently this correspond to the binomial transfer of 
packages of three charges, where the probability of a third 
order process is Pj, = T 3 /(4— 3T) 2 . A similar procedure 
may be applied to the second block Qzb- The result is 
the same. Physically, the two blocks correspond to two 
independent processes which differ by the spin. 

For higher order processes the calculation goes in com- 
plete analogy. The property of a decoupling into two in- 
dependent blocks remains. Further more it is possible to 
the shift the entire x-dependence to the uppermost (or 
the lowest) block. This is achieved by a series of unitary 
operations of the type (1, ..., l,exp(inxT3), 1, ■■■, 1). One 
can easily convince oneself, that for a process of order 
N this gives e. g. the upper-left block K- (N\) and the 
remaining matrix is now independent of x- For example 
a 5th-order process yields 



(50) 



Additionally, the signs of the off-diagonal el- 
ement may be removed by unitary transfor- 
mations. Evaluating the determinant we find 
S( X ) = (2eVt Q /h) In [l + P 5 (e m * - l)] , where 
P 5 = T 5 /(16 - 20T + 5T 2 ) 2 . This expression de- 
scribes binomial transfers of 5 charges with probability 
P 5 - 



f K_(5x) 10 \ 


1 10 


10 1 


10 1 


\ 1 K + (0) J 



Using the above scheme, it is also possible to derive re- 
cursion relations for the probabilities. We find the prob- 
ability for a process of order N 



P 



N 



1 . ((i+^K-i-5)((i-^K-i-5) 



(51) 



The factors a± and 7 are determined from the recursion 
relations 



a$ = l 



T 



T 7„- 



4a£ 



— 7 



*n.-l 

with the initial conditions 



4 OL n _ x a n _^ 



(52) 



7l = VT , a£ = l±?f. (53) 

For general subharmonic voltages 2A/(iV — 1) we find the 
counting statistics 



S( X ) = ^°ln[l + P n (e»*-l) 



where the probabilities are given by 

P 3 = 

P 4 = 

A = 

P 6 = 

P7 = 

Ps = 



(2-T) 2 ' 

(4 - 3T) 2 ' 

(8-8T + T 2 ) 2 ' 

(16-20T + 5T 2 ) 2 ' 

(2-T) 2 (16-16T + T 2 ) 2 ' 

T? 

(64 - 112T + 56T 2 - 7T 3 ) 2 ' 

(T 4 _ 32 T 3 + 160 T 2 _ 256 T + 128) 



(54) 



(55) 



2 ■ 
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Note the limiting cases of these probabilities P/v ~ 
t n i a n-i for T < j and p w = 1 f or J- = 1. 

We can draw several conclusions from the toy model. 
First we obtain simple expressions for the probabilities 
of multiple charge Pn, which are not simple products of 
Andreev reflection probabilities and quasiparticle trans- 
missions, see Eq. (I55|l . Furthermore it is interesting to 
note that by virtue of the unitary transformations we 
can interpret the charge transfer as simultaneous trans- 
mission of N quasiparticles. This explanation does not 
invoke any kind of combined transfer of Cooper pairs and 
quasiparticle. 



B. Full expressions 



Let us now discuss the full expression of the proba- 
bilities P n (E,V) at zero temperature. Since Q has a 
block-tridiagonal form, in order to calculate its determi- 
nant we can use the a recursion technique similar to the 
one describe for the toy model. We define the following 
4x4 matrices 



F± n — Q±n,±n ~ Q±n,±n±2F± n ± 2 Q±n±2,±n ', 71 > 2 

Fq = Qofi — Qo, -2F_ 2 Q-2,0 — Qo, iF 2 Q2,o, (56) 



With these definitions, detQ is simply given by detQ = 
riil-oo detF 2j . In practice, detF„ = 1 if \n\ > A/|eV|. 
This reduces the problem to the calculation of the deter- 
minants of 4 x 4 matrices. 



In the zero-temperature limit one can work out this 
idea analytically, and after very lengthy but straight- 
forward algebra, we obtain the following expressions for 



P^(E,V) = $>_„+ 

i= 

P^E,V) = K 



i=o 



Z' 



1 



n ( r / 4 )i/< 

_k=-n+l+l 
T 



k A \ 2 



Jl ; n>\ 



T 



Y ^-^)-jU±ifBt 2 



Here, we have used again the shorthand g^ ,R (E) = g 



T 



R\2 



(/if) 



A,R 



(E + neV), and defined 



i?<-> A 



(57) 



T 



Z ±n — 1 ± --^-(^±(71+1) ~ 9±n) ~ ~j(f±(n+l)) 5 ±(n+2) > ^ > , 



(58) 



where a = R,A, K = (JI^Li det F-2j)(T[7Li detFV,) and the different functions can be expressed as follows (n > 0) 
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det F± n 



J±n = 



1 ± ^f(.9 a ±n - <£(„-!)) " j(/±„)7^±n , 



n 



Q = A,fl L 



Z2„( 



1±^(*2 



T 



T 



(JJ det F ±(n+2j .)) <^ — (<^„ - 3 £„ 

3=1 I 



(59) 



a \2 



(ff±n " 5±0-l))) _ 'j(/±n) 



r 



T, 



ryK ryA I J^l 



A |2 



T ~J\f±n ~ f±n) [f±n,Z±n + f±n^±n\ 



Notice that, since at zero temperature the charge only 
flows in one direction, only the P n with n > survive. 
It is worth stressing that the full information about the 
transport properties of superconducting point contacts is 
encoded in these probabilities. Let us also remark that 
P n {E, V) are positive numbers bounded between and 1, 
and fulfill the normalization condition J2 n Pn(E, V) = 1. 
Thus, we see that for the finite bias dc transport, where 
the superconducting phase does not play any role, there is 
no problem with the typical interpretation of P n as prob- 
abilities [I3. Moreover, although at a first glance the 
expressions of Eqs. (|57I59|I look complicated, they can 
be easily computed and provide the most efficient way 
to calculate the transport properties of these contacts. 
In practice, to determine the functions B A ' R and det F n , 
one can use the boundary condition B A ' R — det F n = 1 
for |n| >A/|eV|. 

In view of Eqs. (|57I59|I the probabilities P n can be in- 
terpreted in the following way. P n is the probability of 
a MAR of order n, where a quasiparticle in an occupied 
state at energy E is transmitted to an empty state at 
energy E + neV. The typical structure of the expres- 
sion for this probability consists of the product of three 
terms. First, Jo gives the probability to inject the incom- 
ing quasiparticle at energy E. The term n£=i (T/^)\fk\ 2 
describes the cascade of n — 1 Andreev reflections, in 
which an electron is reflected as a hole and vice versa, 
gaining an energy eV in each reflection. Finally, J n gives 
the probability to inject a quasiparticle in an empty state 
at energy E + neV. This interpretation is illustrated in 
Fig.^ where we show the first four processes for BCS su- 
perconductors. The product of the determinants in the 
expression of J„ (see Eq. |JS5J) describes the possibility 
that a quasiparticle makes an excursion to energies below 
E or above E + neV. In the tunnel regime this possibility 
is very unlikely and at perfect transparency is forbidden. 
For this reason the expressions of the MAR probabilities 
simplify a lot in these two limits, as we discuss in the 
next paragraphs. 

In the tunnel regime a perturbative calculation yields 
(n > 1) 



P»(r«i) 



rpn 



where p(E) is the reservoir density of states. If we use 
this result in the current expression (see below), we re- 
cover exactly the result of the multiparticle tunneling 
theory of Schrieffer and Wilkins 9]. As we mentioned 
in the introduction, the expression above leads to diver- 
gences in the current, which shows that this problem 
is non-perturbative in the transmission. Thus, even at 
low transparencies one has to use the full expression of 
Eqs. (|57I59() . where the mentioned divergences are regu- 
larized in a natural manner. 

For perfect transparency (T = 1), the absence of nor- 
mal backscattering makes the expressions of the proba- 
bilities P n (E, V) much simpler, and one can show that 
they can be written as (n > 1) 



n-1 

P n (T = 1) = £(Ha- 



n+l\ ) 



1-1 

n 

_k=-n+l+l 



Ofc 



-POP, 



*:=i 



t\ 2 



(60) 



(i-hl 2 ), 

(61) 

where a(E) is the Andreev reflection coefficient defined 
as a(E) = -if R (E)/ [l + g R {E)} , and a n = a(E + neV). 
As can be seen in Eq. (|61|) , a quasiparticle can only move 
upwards in energy due to the absence of normal reflec- 
tion. If we use this expression in the current formula 
we recover the result obtain by Klapwijk, Blonder and 
Tinkham H3 for T = 1. 



IV. APPLICATION TO DIFFERENT 
SITUATIONS 

As explained in the previous section, with the expres- 
sion of the MAR probabilities we can easily describe 
many different transport properties. Moreover, notice 
that so far we have not made any assumption about the 
leads Green's functions g A ' R and f A:R entering in the 
expressions of P n {E, V). Therefore, these expressions al- 
low us to address a great variety of situations. In this 
section we analyze the zero-temperature transport prop- 
erties of three different situations: (i) a contact between 
BCS superconductors, (ii) a contact between supercon- 
ductor under the influence of pair-breaking mechanisms 
and (iii) a short diffusive SNS contact, where N is a nor- 
mal disordered region shorter than the superconducting 
coherence length. 
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FIG. 2: Current, shot noise and third cumulant at zero tem- 
perature as a function of the voltage for BCS superconduc- 
tors of gap A. The different curves correspond to different 
transmission coefficients as indicated in the panels. Here, 
Gn — (2e 2 /h)T is the normal state conductance. 



A. BCS superconductors 

Let us start analyzing the most standard situation, 
namely a contact between two BCS superconductors with 
a gap A. In this case f A > R = iA/y/(E =F iv) 2 ~ A 2 , 
where r\ = + , and g A ' R follows from normalization. 
As mentioned in the introduction the current and noise 
of such a contact have been thoroughly studied both 
theoretic ally ITM Flrl JTA E IH H^j and experimen- 
tally p^EllSnmiLIS^. Our goal here is to show how the 
knowledge of the FCS provides a new and deeper insight 
into the different transport properties. 

In Fig. [2 we show the first three cumulants of the 
charge transfer distribution: current, shot noise and 
skewness (third cumulant). Let us discuss their most 
remarkable features, (i) The current exhibits the so- 
called subharmonic gap structure, as discussed in the 
introduction. This subgap structure evolves from a step- 
like behavior for low transmission to its disappearance 




FIG. 3: (a) Second cumulant and (b) third cumulant at zero 
temperature for BCS superconductors. Both are normalized 
to the first cumulant (the average current). The transmissions 
are indicated in the plots. 



at perfect transparency, (ii) The shot noise in the sub- 
gap region can be much larger than the Poisson noise 
{Si^poisson = 2eJ). Moreover, in the tunneling regime 
the effective charge defined as the ratio q = Sj/27 is 
quantized in units of the electron charge: q(V)/e — 
1 + Int(2A/eV). This is illustrated in Fig. where 
the ratios C%IC\ and C3/C1 are shown as a function 
of the voltage, (iii) As shown in Fig. |3 the third cu- 
mulant at low transmissions is described by C3 = q 2 Ci, 
where again q is the quantized effective charge defined 
above. For higher transmissions this cumulant is neg- 
ative at high voltage as in the normal state, where 
C* 3 = (t /h)T(l - T)(l - 2T)eV, but it becomes pos- 
itive at low bias, and after this sign change there is a 
huge increase of the ratio C3/C1. 

The features described in the previous paragraph can 
be easily understood with the help of an analysis of the 
probabilities P n (E,V). To give an idea about them, in 
Fig. 2] we have plotted their average value defined as 
P n (V) = (1/eV) f° V dE P„(E,V) for two very differ- 
ent transmissions. First of all, notice that, no matter 
what the transmission is, the probability of an n-order 
MAR has a threshold voltage eV n — 2A/n, below which 
the process is forbidden. When V > V n an n-order MAR 
gives a new contribution to the transport, which is fi- 
nally the explanation of the subharmonic gap structure. 
On the other hand, the big difference between the tun- 
neling regime and perfect transparency can be explained 
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FIG. 4: Average MAR probabilities P n (V) = 

(1/eV) J_ e dE P„(E,V) as a function of voltage for a 
contact between BCS superconductors at zero temperature. 
The two panels correspond to two different transmissions. 
The index of the processes is indicated in the plots. Notice 
the logarithmic scale in the panel (a). 



as follows. At low transparency there are two factors 
that make the subgap structure so pronounced. First, at 
V n the n-order MAR is a process that connects the two 
gap edges, where the BCS density of states diverges (see 
Eq. H60|l). This fact, together of course with its higher 
probability, implies that this MAR rapidily dominates 
the shape of the I-V curves giving rise to a non-linearity 
at V n . Second, at V n there is a huge enhancement of 
the probabilities of the MARs of order m > n. This is 
due to the fact that precisely at V„ the MAR trajecto- 
ries can connect both gap edges, which as can be seen 
in Eq. (|60|l increases enormously their probability. At 
perfect transparency, the MAR probabilities do not ex- 
hibit any abrupt feature (see Fig. 0Jd) . This is due to the 
fact that the BCS density of states is renormalized, and 
in particular, the divergences disappear (see Eq. H61[) ). 
This fact explains naturally why the subharmonic gap 
structure is completely washed out at T = 1. 

Another interesting feature of the MAR probabilities 
occurs at low transparencies. As one can see in Fig. |3Ji, 
at a voltage 2A/n < eV < 2A/(n - 1) the MAR of order 
n has a much higher probability than the other MARs. 
This means that in this voltage window the n-order MAR 
clearly dominates the transport properties and the charge 
is predominantly transferred in packets of ne. This fact 
explains the charge quantization in the tunnel regime ob- 
served both in Ci and C3 (see Fig- EJ- More generally, 
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this fact implies that at low transparencies the multi- 
nomial distribution of Ea. (|19J) becomes Poissonian, and 
in this limit all the cumulants are proportional to the 
current: C n — (q(V)/e) n Ci, where q(V) is the voltage- 
dependent quantized charge. When the transmission is 
not very low, there are always several MARs that give a 
significant contribution to the transport at every voltage 
(see Fig.0J>). This explains why the charge is in general 
not quantized. 

The explanation for the sign change of C3 at low bias 
and high transparencies can be found in Eq. (23). In 
order to get a positive value for C3, one needs the first 
two terms in Eq. (23) to dominate, which happens when 
P n <C 1. This is precisely what happens at low bias, 
where the MAR probabilities are rather small. On the 
other hand, the huge enhancement after the sign change 
is due to fact that n, the charge transferred by these 
MARs, is indeed huge at low bias. 

Finally, at T = 1 the cumulants C n (with n > 1) do 
not completely vanish due to the fact that at a given 
voltage different MARs give a significant contribution, 
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FIG. 6: Zero temperature current-voltage characteristics for 
superconductors with a depairing energy F in units of Ao. 
The current and the voltage have been normalized with the 
order parameter A at the corresponding F. The different 
panels correspond to different transmissions values. 



and therefore their probability is smaller than one (see 
Fig. lib)). 



B. Pair-breaking mechanisms 

It is well-known that there are many mechanisms that 
can lead to pair-breaking effects, which modify the quasi- 
particle spectrum of a superconductor. Typical examples 
are a magnetic field, supercurrents or magnetic impuri- 
ties. It was shown in the 1960's that for diffusive su- 
perconductors various pair-breaking mechanisms can be 
described in a unified manner with a single parameter 
F, the depairing energy, which describes the strength 
of the pair- breaking J65| |. The only difference between 
these mechanisms is contained in the microscopic ex- 
pression of r. For instance, for thin a film of thick- 
ness d much smaller than the superconducting coher- 
ence length in a magnetic field H parallel to the film 
T = De 2 d 2 H 2 /(6hc 2 ), where D is the diffusion con- 
stant. In these situations the energy-dependent retarded 
Green's function can be calculated from 65] 




FIG. 7: Current contribution of processes n = 1,2,... for 
T — 1 as a function of voltage for superconductors with a 
depairing energy F in units of Ao. The current and the volt- 
age have been normalized with the order parameter A at the 
corresponding F. The order of the processes is indicated in 
the plots. 




FIG. 8: Zero temperature noise for superconductors with 
a depairing energy F in units of Ao. The current and the 
voltage have been normalized with the order parameter A 
at the corresponding T. The different panels correspond to 
different transmissions values. 
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(62) 



Here, A is the order parameter, which is in this case dif- 
fers from the spectral gap and it has to be determined 
self-consistently [66|. For small T the pair-breaking 
mechanisms result in a smearing of the BCS singular- 
ities in the density of states and in a suppression of 
the spectral energy gap A g to a reduced value A g = 

A [l — (r/A) 2 / 3 ] . The gap disappears completely at 
r ~ 0.45Ao, where Ao is the order parameter in the 
absence of pair-breaking. The gapless superconductivity 



survives until the critical value r^ = 0.5Ao. This behav- 
ior is illustrated in Fig.^a), where we show the density 
of states as a function of energy for different values of T 
in units of Ao. In Fig.[S{b) one can see the evolution of 
the order parameter and spectral gap with the depairing 
energy. 

Let us discuss now how this modified density of states 
is reflected in the transport properties. In Fig. we 
show I-Vs for different transmissions and different val- 
ues of the depairing energy. The most noticeable fea- 
tures are: (i) the subharmonic gap structure is shifted 
to voltages eV — 2A g /n, and (ii) the subgap structure 
progressively disappears as the pair-breaking strength is 
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FIG. 9: Zero temperature third cumulant for superconduc- 
tors with a depairing energy F in units of Ao . The current and 
the voltage have been normalized with the order parameter 
A at the corresponding F. The different panels correspond to 
different transmissions values. 



increased. These features are simple consequences of 
the evolution of the density of states with T. Anyway, 
one can get a further insight by analyzing the contri- 
bution to the current of the individual MAR processes: 
I n = (2e/ft.) J Q e dE P n (E, V). These quantities are plot- 
ted in Fig. for T = 1. As one can see, the threshold 
voltage for a n-order MAR is now eV n = A g /2n as a 
consequence of the reduced spectral gap. As the gap di- 
minishes, the processes of lowest order dominate the I-Vs 
even at low bias. It is interesting to notice that even in a 
gapless situation (r = 0.475) there is a finite contribution 
of the MARs. It is worth mentioning that in Refs. |67J 
and [68j the type of theory described here accounted for 
the magnetic field dependence of the I-Vs of atomic con- 
tacts. 

Let us turn now our attention to the second and third 
cumulants that can be seen in Fig. |S] and Fig. El respec- 
tively. As in the case of the current, the subharmonic gap 
structure is shifted and smoothed as the gap evolves with 
r. Moreover, one can notice that for high transparencies 
and in the subgap region there is a great reduction of 
both cumulants as T increases. This is a consequence of 
the fact that low order MARs dominate even at low bias, 
which in practice means that the charge transferred at 
these voltages is on average not very big. 



C. Diffusive SNS contacts 

So far we have discussed the case of a single channel 
contact. The results are trivially generalized to the mul- 
tichannel case by introducing a sum over the conduction 
channels. In this subsection we briefly address the case 
of a short diffusive SNS junction with a large number 
of transmission channels and diffusive electron transport 
in the normal N region. The superconducting leads are 
considered as BCS superconductors. In this case, the dis- 
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FIG. 10: Zero temperature transport properties of a short 
diffusive SNS junction, (a) First three cumulants: current in 
units of (Gn A/e), shot noise in units of (2Gn A) and the third 
cumulant in units of (GjvAto/ftGo). (b) Current contribution 
of the different processes, (c) Ratio C2/C1 as a function of 
the inverse of the voltage, (d) Ratio C3/C1 as as a function 
of the inverse of the voltage. 



tribution of transmission coefficients is continuous, and 
it is characterized by the density function p(T), which 
has the well-known bimodal form 691 



P(T) = 



G 



y 



1 



2GoTVT^T 



(63) 



where Gn is the normal-state conductance of the N re- 
gion ad Go = 2e 2 /h is the conductance quantum. Then, 
the different cumulants can be calculated from the single- 
channel results C n (T) as follows 



L „ 



dTp(T)C n (T). 



(64) 



In Fig. llOf a^) we show the first three cumulants for this 
SNS system. Both the current and the noise have pre- 
viously discussed in the literature [23, IZfl > an( i nere we 
recover these results. Both quantities exhibit a subhar- 
monic gap structure which is a result of the competition 
of channels with different transparencies. Again, this 
structure can be understood by analyzing the individ- 
ual contributions to the current of the different MARs, 
see Fig. llOf b). As one can see, at every voltage there 
are several processes giving a significant contributions, 
which makes that subgap structure much smoother than 
in the single-channel case. This fact also explains the 
absence of the charge quantization in this multichannel 
case. This is illustrated in Fig. llUf c). where we show 
the ratio C2/C1 as a measure of the effective charge. 
Notice that at low bias this effective charge grows as 
(1/V) as obtained in Ref. |23. In this regime the nu- 
merical results can be approximately described by the 
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FIG. 11: Current-voltage characteristics at finite tempera- 
ture for BCS superconductors. The temperature is in units of 
the critical temperature Tc- The current and the voltage are 
normalized with the temperature-dependent gap. The differ- 
ent panels correspond to different transmission values. 



following linear function: C 2 /Ci = 0.31(2A/eV) + 0.55. 
On the other hand, the third cumulant exhibits a huge 
increase at low voltages |35j. In particular, as shown 
in Fig. HUr d). the ratio C3/C1 grows as (l/V) 2 at low 
bias. In this regime the ratio can be approximated by 
C3/C1 = 0.05(2A/eV) 2 + 0.5. 



V. TRANSPORT PROPERTIES AT FINITE 
TEMPERATURES 

So far we have discussed the transport properties of 
superconducting point contacts at zero temperature. In 
this section we shall investigate the role of the tempera- 
ture, which we shall denote as T e . We focus our attention 
to the case of a single channel contact between BCS su- 
perconductors. At finite temperature it is not easy to 
determine analytically the probabilities P n (E,V), and 
in this case we have calculated them numerically. The 
idea goes as follows. According to Eq. (18) we need to 
calculate the coefficients Pn(E, V), which are simply the 
Fourier coefficients of the series in Eq. (18), i.e. 



K(e, v ) = ^J* d x e ~ mx det Q(x) 



(65) 



Finally, detQ(x) is calculated numerically. Of course, if 
one is only interested in the different cumulant, one can 
easily calculate them by taking the numerical derivative 
oftheCG F see Eq.El 

In Figs. lllfT^l and 1131 we show the current, noise and 
third cumulant, respectively, for different transmission 
and temperatures ranging from zero to the critical one. 
Notice that in order to get rid of the trivial temperature 
dependence due to the decrease of the gap we have nor- 
malized the voltage by the temperature-dependent gap 
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FIG. 12: Finite temperature noise for BCS superconductors. 
The temperature is normalized with the critical temperature 
Tc- The different panels correspond to different transmis- 
sion values. The voltage is normalized with the temperature- 
dependent gap, and the current with the zero-temperature 
gap. Note that the scaling is different from the other plots in 
this section. 
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FIG. 13: Finite temperature thrid cumulant for BCS super- 
conductors. The temperature is normalized with the critical 
temperature Tc- The different panels correspond to different 
transmission values. The third cumulant and the voltage are 
normalized with the temperature-dependent gap. 



A(T e ). As it can be seen in Fig. II II the temperature pro- 
gressively smoothes the SGS and increases the current for 
low transmissions. These are simple consequences of the 
thermal excitation of quasiparticles. For higher trans- 
missions the temperature has the opposite effect (see 
the lower two panels in Fig. Illjl . The current decreases 
with increasing temperature and approaches the normal 
state current-voltage characteristic from above. At the 
same time the excess current, i. e. I(V » A/e) — GnV, 
vanishes obviously. So in short, by increasing the tem- 
perature high-order Andreev reflections contribute less 
to the current, which is dominated by thermally acti- 
vated direct quasiparticle tunneling. This behavior is 
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FIG. 14: Average MAR probabilities P n (V) = 

(1/eV) L e d-E P n {E, V) at finite temperature as a function 
of voltage for a contact between BCS superconductors with 
transmission T — 0.95. The four panels correspond to dif- 
ferent temperatures T e expressed in units of the critical tem- 
perature Tc- The index of the processes is indicated in the 
plots. 



clearly illustrated in Fig. 1141 where we show the evo- 
lution with the temperature of the average probability 
P n (V) = (1/eV) f* V dE P n (E,V) of different processes 
for a contact with transmission T = 0.95. Notice that 
we only show the first electron processes that give a pos- 
itive contribution to the current. Remember that at fi- 
nite temperature there are also hole processes that give 
a negative contribution to the current, the magnitude of 
which is still much smaller than the one of the electron 
processes in the shot noise limit eV ^S> ksT. At van- 
ishing voltages, of course, P n — P- n as required by the 
fluctuation-dissipation theorem. 

In Fig.^Jone can observe the following important fea- 
tures. First, at finite temperature the different processes 
do not have any finite threshold voltage, and they can 
contribute down to zero bias due to thermal activation. 
Second, as the temperature increases the probability of 
the single quasiparticle processes is greatly enhanced in- 
side the gap. This fact results in a reduction of the av- 
erage effective charge transmitted through the contact. 
Finally, notice that although the MAR probabilities are 
reduced inside the gap at finite temperature, high-order 
processes can give a significant contribution to the trans- 
port even at voltages larger than the gap at the corre- 
sponding temperature. This is clearly at variance with 
the zero temperature case. To understand this behav- 
ior, let us recall that the total voltage gain for an order 
n process is neV, which means essentially that higher 
order processes can start well below the gap and end 
well above the gap. Now, at finite temperature e.g. the 
end states above the gap are filled with finite probabil- 



ity f(E + neV), assuming that the process has started 
with a quasiparticle at energy E. A certain process can 
only happen if its final state is empty. This gives a factor 
1 — f(E + neV), which enhances the chance for higher 
order processes, since they have to end up at higher ener- 
gies, for which this factor is larger . On the other hand, 
a similar argument can be made about the initial state, 
which has to be filled for the process to take place. Again, 
this is more likely for higher order processes, since they 
can emerge from energies well below the gap, which are 
completely filled also at finite temperature. 

It is interesting to discuss the qualitative different tem- 
perature behavior of the second and third cumulants. 
The noise exhibits a transition from pure shot noise at 
zero temperature to thermal noise when the tempera- 
ture is larger than the voltage. As it can be seen in 
Fig. ^1 this transition is reflected in a saturation of the 
noise at low bias to a finite value, which is given by the 
fluctuation-dissipation theorem. It is interesting to note, 
that the noise decreases as a function of voltage in the 
transition region from thermal to shot noise also for rel- 
atively small transmissions. Such a behavior can be at- 
tributed to the multinomial distribution. Interestingly, 
from Eq. (|28l) we see that the correlations between pro- 
cesses of orders with opposite sign (e. g m = —n) tend 
to increase the noise. As these terms appear only if the 
respective probabilities are non-negligible, the reduction 
of noise below the thermal level can be understood as 
consequence of the vanishing cross correlations between 
processes of orders with different signs. 

The temperature dependence of the third cumulant is 
very interesting. First we recall that the third cumulant 
vanishes at zero voltage for any temperature (as all odd 
cumulants do). In Ref. J4J| the temperature dependence 
of the third cumulant for a quantum contact between 
normal metals was calculated. It was shown, that an in- 
creasing transparency has quite a dramatic effect on the 
third cumulant. For a tunnel junction (i. e. for small 
transmission) C3 is independent of the temperature and 
it is simply equal to the q 2 C\. However, this is interest- 
ing because it allows a direct measurement of the charge 
q transfered in an elementary event even for voltages be- 
low the shot noise limit. Note, that this relation holds 
also for non-linear current-voltage characteristics, since 
it is a consequence of the bidirectional Poisson distribu- 
tion in this limit. The effects of a finite transparency are 
even more dramatic. The third cumulant has a marked 
temperature dependence, crossing over from a FI depen- 
dence, where F = 1 — T is the Fano factor, to a novel 
high-temperature dependence ~ FI(1 — 2T), which can 
even become negative for T > 1/2. In view of these find- 
ings, we will now discuss our results for the temperature 
dependence of the third cumulant of a superconducting 
point contact. 

First, we note that in Fig. ^] C3 has a temperature 
dependence even in the tunnel regime. As explained in 
the previous paragraph, this in contrast with the normal 
state, where C3 is almost independent of the tempera- 
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FIG. 15: Ratio C3/C1 for two different temperatures as a 
function of the voltage for a contact between BCS supercon- 
ductors. The different curves correspond to different trans- 
missions as indicated in the plots. 



ture, as it has been discussed theoretically in Ref. J45J 
and observed experimentally in Ref. |71j . In our case the 
temperature dependence is due to the change in the MAR 
probabilities caused by the thermal activation. As ex- 
plained above, the thermal activation enhances the prob- 
ability of the tunneling of single quasiparticles inside the 
gap, which in turn reduces the average effective charge. 
A consequence of this fact is the great reduction of the 
ratio C3/C1 as the temperature increases. This is illus- 
trated in Fig 1151 This reduction is specially dramatic in 
the subgap region for high transparencies, as it can be 
seen directly in Fig. ^J 



VI. CONCLUSIONS 

We have presented a detailed analysis of the full count- 
ing statistics in superconducting point contacts at finite 
bias voltage. We have demonstrated that the charge 
transfer in these systems is described by a multinomial 
distribution of processes, in which multiple charges ne 
(with n = 1, 2, 3, ..., 20, ...) are transferred through the 
contact. These processes are nothing but multiple An- 
dreev reflections. The knowlegde of the full counting 
statistics allows us to obtain the probabilities of the in- 



dividual MARs, providing so a deep insight into the elec- 
tronic transport of these junctions. From the knowledge 
of these probabilities one can easily calculate not only the 
current or the noise, but all the cumulants of the current 
distribution. We have also shown that one can obtain 
analytical expressions for the MAR probabilities at zero 
temperature, which provides the most efficient method 
to calculate the transport properties of these contacts. 
Moreover, the FCS approach allows us to describe a great 
variety of situations in a unified manner. 

In this sense, we have addressed different situations 
such as contacts between BCS superconductors, junc- 
tions between superconductors where a pair-breaking 
mechanism ist acting or short diffusive SNS contacts. 
We have also discussed the temperature dependence of 
the first cumulants and illustrated their peculiarities as 
compared with the normal case. It is also worth men- 
tioning that the formalism developed in this work can 
be easily applied to other situations not addressed here 
such as point contacts with proximity-effect supercon- 
ductors [72| and Josephson junctions of unconventional 
superconductors |73ll74|. 

From the full counting statistics view, we have found a 
new distribution occuring in superconducting point con- 
tacts. The statistics takes the form of a multinomial 
distribution of charge transfers of all orders, which are 
allowed by the applied bias voltage. We have shown, 
that the limit of opaque contacts provides an interesting 
situation, in which Poissonian statistics makes it possi- 
ble to observe multiple charge transfers in a direct man- 
ner. Furthermore, we have discussed consequences of the 
multinomial statistics of charge transfers of different sizes 
at the same time. For example, an open contact has a 
finite noise due to the presence of different MAR pro- 
cesses at the same time. The temperature dependence 
of the counting statistics provides a new insight in the 
transport characteristic, since we have shown that higher 
order Andreev processes contribute also at voltages much 
larger than the superconducting gap. 

Finally we remark, that the FCS approach provides a 
fresh view of the electronic transport of superconduct- 
ing point contacts and it is seems to be a natural choice 
for future theoretical analyses. On the other side, super- 
conducting contacts show an interesting new counting 
statistics, namely a multinomial distribution, and we ex- 
pect further intersting results in other superconducting 
systems out of equilibrium. 
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